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Thermal energy storage is needed to improve the efficiency of solar thermal energy applications (STEA) 
and to eliminate the mismatch between energy supply and energy demand. Among the thermal energy 
storages, the latent heat thermal energy storage (LHTES) has gained much attention because of its high- 
energy densities per unit mass/volume at nearly constant temperatures. This review presents previous 
studies on the numerical modeling of phase change materials (PCMs) through a commercial computa- 
tional fluid dynamic (CFD) software and self-developed programming to study the heat transfer 
phenomena in PCMs. The CFD (Fluent) software is successively used to simulate the application of 


Keywords: 
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CFD PCMs in different engineering applications, including electronic cooling technology, building thermal 
puent storage, and heating, ventilation, air conditioning (HVAC). Using CFD software to design LHTES is 
believed to be an effective way to save money and time and to deliver optimization tools for maximum 
efficiency of STEAs. 

© 2012 Elsevier Ltd. All rights reserved. 
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1. Introduction ranges. PCM thermal storage plays a key role in improving energy 


efficiency. It limits the discrepancy between the energy supply 


Researchers intensively studied the thermal energy storage 
of PCMs for the last three decades because of the latter’s high 
thermal energy densities per unit volume/mass and their applic- 
ability to different engineering fields using wide temperature 
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and the energy demand of solar thermal energy applications 
(STEAs), particularly when the STEA operation strategy depends 
solely on solar energy as a main source. PCM thermal storage 
indicates high performance and dependability with the advan- 
tages of high storage capacity and nearly constant thermal energy 
[1]. Most STEAs need constant or near-constant temperature for 
high-efficiency strategies. Using latent heat thermal energy sto- 
rage (LHTES) as thermal energy storage can provide the required 
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Nomenclature 

C Mushy zone constant (kg/m? s) 
Cp Specific heat of PCM, (J/kg) 

g gravity acceleration, (m/s?) 

h sensible enthalpy(J/kg) 

H enthalpy (J/kg) 

L latent heat fusion (J/kg) 

k thermal conductivity (W/m K) 
T Temperature (°C or K) 

u velocity (m/s) 

Greek letters 


P fluid density (kg/m?) 

y liquid fraction 

B volumetric expansion coefficient (1/K) 
H Dynamic viscosity (kg/m s) 

& Phase volume fraction 

Subscripts 

ref Reference 

S Solidus of the phase change materials 
l Liquidus of the phase change material 


constant temperature that matches the melting temperature of 
the PCMs. PCMs are used in different engineering fields, such as in 
the following: the thermal storage of building structures [3-7]; 
building equipment, including domestic hot water; heating and 
cooling systems [8,9]; electronic products [10-12]; drying tech- 
nology [13], waste heat recovery [14]; refrigeration and cold 
storage [15-18]; solar air collectors[19]; and solar cookers [20]. 
Using a CFD software to design a LHTES is expected to be an 
effective way to save money and time and to deliver optimization 
tools for maximum efficiency of STEAs. 


2. Numerical solution of PCMs 


The mathematical formulation of a phase transient known 
as a phase change or moving boundary is governed by a partial 
deferential equation that can be solved either analytically or 
numerically. The analytical solution of PCMs is problematic 
because of the nonlinear phase front interfaces, complex geome- 
tries, and nonstandard boundary condition; the few analytical 
studies available are on 1D cases with regular geometries and 
a standard boundary condition. Ozisik [21] reported that the 
numerical method for solving PCMs can be categorized as fixed- 
grid, variable grid, front-fixing, adaptive grid generation, and 
enthalpy methods. 

Predicting the behavior of phase change systems is difficult 
because of its inherent non-linear nature at moving interfaces, for 
which the displacement rate is controlled by latent heat lost or 
absorbed at the boundary [22]. The heat transfer phenomena in 
solid-liquid PCMs can be analyzed using two main methods: the 
temperature-based and enthalpy-based methods. In the first 
method, temperature is considered a sole dependent variable. 
The energy conservation equations for the solid and liquid are 
written separately; thus, the solid-liquid interface position can be 
tracked explicitly to achieve an accurate solution for the pro- 
blems. 


Ts, _ 
als k= 


oTi 
n ay k,+ pLkvy, (1) 


where T, denotes the temperature in the solid phase, T, denotes 
the temperature in the liquid phase, k, is the thermal conductivity 
of the solid phase, kı is the thermal conductivity of the liquid 
phase, n is the unit normal vector to the interface, and v, is the 
normal component of the velocity of the interface. L is the latent 
heat of freezing, as shown in Fig. 1. 

In the second method, the solid-liquid interface position need 
not be tracked. Researchers often use the enthalpy formulation 
because of the following advantages: (1) the governing equations 
are similar to the single-phase Eq.; (2) no explicit conditions need 


to be satisfied at the solid-liquid interface; (3) the enthalpy 
formulation involves the solution within a mushy zone, involving 
both solid and liquid materials, between the two standard phases; 
and (4) the phase change problem can be solved more easily [23]. 
APD LY. (pvH) = V-(kAT) +S (2) 

Where T denotes the temperature, k is the thermal conductivity, 
p is the density of the PCM, 7 is the fluid velocity, and H is the 
enthalpy, S is the source term. 

Dutil et al. [22] presented an intensive mathematical and 
numerical review of the PCM application based on the first and 
second laws of thermodynamics. Using 252 references, they 
determined the mathematical and numerical methods applied 
to solve heat transfer problems involving PCMs for thermal 
energy storage, the mathematical fundamentals of PCMs, and 
the different application geometries and applications. Verma et al. 
[24] also introduced other PCM mathematical reviews. 

Recently, researchers used the Fluent software by ANSYS to 
simulate melting and solidification in engineering problems. 
Other software that can be used to simulate the PCM process 
include COMSOL Multiphysics and Star-CMM +. However, Fluent 
is preferred by most researchers for melting and solidifying PCMs. 
Conversely, some researchers self-developed a program using 
computational language (C+4, Fortran, Matlab) to study the heat 
transfer phenomena in PCMs. Table 1 summarizes some of the 
self-developed programs for different PCM geometries. 


2.1. Fluent program 


The Fluent software by ANSYS is a computational fluid 
dynamic (CFD) program used successfully to simulate different 
engineering problems. This software has a specific model that can 
simulate a range of different melting and solidification problems 
in engineering, including casting, melting, crystal growth, and 


a 


Fig. 1. Solid-liquid interface for a multidimensional situation [21]. 


Table 1 


Self-developed numerical model for different PCM geometries. 
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Ref PCM Geometry Dimension Discretization Numerical solution methods Validation 
methods 
[66] Gallium Rectangular 2D (X, Y) FDM Enthalpy method-based control volume approach; the power law used for Experimental 
discretization and semi-implicit method for pressure-linked equations (SIMPLE) was 
adopted for the pressure-velocity couple and for the tridiagonal matrix algorithm 
(TDMA) solver to solve the algebraic equations 
[67] n-Octadecane Rectangular 1D (X)2D FDM An enthalpy formation-based control volume approach, the fully implicit method, and [68,69] 
(X, Y) TDMA were used to solve the algebraic equations 
[70] Paraffin wax PCM bags 3D (X,Y,Z) FDM The interpolating cubic spline function method was used to determine the specific [71] 
RII-56 In duct heat in each time step of the temperature calculations; the control volume approach 
and the backward Euler scheme were chosen for the temperature discretization in 
time 
[72] Paraffin RT30 Cylindrical 2D (R, X) FVM The power law used for discretization and SIMPLE were adopted for pressure-velocity Experimental 
coupling 
[73] n-Octadecane Cylindrical 2D (R, X) FDM Enthalpy methods and the TDMA solver were used to solve the equations Experimental 
[74] 4 PCM Cylindrical 2D (R, X) FDM Fixed control volume and enthalpy methods - 
[75] RT60 Cylindrical 2D (R, ©) FDM A control volume-based implicit form was used. All equations were solved Experimental 
simultaneously using the Gauss-Seidel iterative method 
[76] Paraffin 130/ Cylindrical 2D (R, ©) FDM The enthalpy formulation approach, the control volume method, and the alternating Experimental 
135 Type 1 direction scheme were used to discretize the basic equations 
[77] Paraffin wax Cylindrical 2D (R, ©) FDM The control volume approach, the power law used for discretization, and the SIMPLE [78,79] 
algorithm were adopted for pressure-velocity coupling 
[80] n-Hexacosane Cylindrical 1D (R) - A temperature and thermal resistance iteration method was developed to analyze Experimental 
PCM solidification to obtain the solutions 
[81] 4PCM Cylindrical 2D (R, X) FDM Enthalpy methods were used; the equations were solved simultaneously using the 
Gauss-Seidel iterative method 
[82] 2PCM Cylindrical 2D(R,X) FDM The enthalpy method and the TDMA were used to solve the resulting algebraic [73] 
equations solved using an iteration procedure 
[83] 5 PCM Cylindrical 2D (R, X) FEM Standard Galerkin finite element method was used - 
[84] 3 PCM Cylindrical 2D (R, X) FDM Using the enthalpy method, the equations were solved by adopting the alternating [85] 
direction implicit (ADI) method and the TDMA 
[86] Water Spherical 1D (R) FDM An implicit FDM approach and a moving-grid scheme were used. The problem domain Experimental 
was divided into two regions, solid and liquid, both of which were separated by the 
interface. All algebraic equations were solved simultaneously 
[87] Water, water Spherical 1D (R) FDM The heat delivered by the capsule could also be calculated by the thermal resistance Experimental 
Glycol circuit. The resulting equation and associated boundary conditions were treated using 
mixtures a moving grid method 


FDM; Finite difference method, FEM; Finite elements method, FVM; Finite volume method. 


solidification. The program can be used to solve the phase change 
that occurs at a single temperature (pure metals) or over a range 
of temperatures (mixture, alloy, and so on). The applications and 
limitations of Fluent can be found in Ref. [25]. 

To begin the Fluent software simulation, the physical engi- 
neering problem is drawn and meshed in a specific geometric 
modeling using mesh generation tools that include the Fluent 
software (Gambit, workbench). The mesh generation tools can be 
imported from other programs like Auto CAD. After the physical 
configuration is drawn and meshed, the boundary layers and 
zone types are defined, and the mesh is exported to the Fluent 
software. 

Different grid sizes and time steps should be applied to the 
numerical model to ensure that the numerical results are inde- 
pendent of the parameters. Small grid sizes and time steps are 
preferred for a short simulation time in the computer. 


2.2. Mathematical formulation for Fluent 


The mathematical equations used to solve the solidification 
and melting models in Fluent depend on the enthalpy-—porosity 
technique [26-28] and on the finite volume methods. In the 
former, the melt interface is not tracked explicitly. A quantity 
called liquid fraction, which indicates the fraction of the cell 
volume in liquid form, is associated with each cell in the domain. 
The liquid fraction is computed at each iteration based on 
enthalpy balance. The mushy zone is a region wherein the liquid 
fraction lies between O and 1. The mushy zone is modeled as a 
“pseudo” porous medium in which the porosity decreases from 


1 to O as the material solidifies. When the material has fully 
solidified in a cell, the porosity becomes zero, resulting in the 
drop of velocities to zero. In this section, an overview of the 
solidification/melting theory is given. For details on the enthalpy- 
porosity method, refer to Voller and Prakash [26]. 

The energy equation is expressed as 


6(PH) + V (PTH) = V(kVT)+5S (3) 


where p is the density of the PCM, 7 is the fluid velocity, k is the 
thermal conductivity, H is the enthalpy, and S is the source term. 
The sensible enthalpy can be expressed as 


T 
h=ħet | cpAT, (4) 
ref 
and H can be defined as 
H=h+AH, (5) 


where Iyer is the reference enthalpy at the reference temperature 
Tref, cp is the specific heat, AH is the latent heat content that may 
change between zero (solid) and 1 (liquid), L is the latent heat of 
the PCM, and y is the liquid fraction that occurs during the phase 
change between the solid and liquid state when the temperature 
is T; >T > Ts. Thus, y may be written as 


y=AH/L, (6) 
0 if T<T,; 

y= 1 if T>T, (7) 
(T—Ts)/(T\-Ts) if Tı >T> Ts 
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The source term S; in momentum equation, is defined as 
v 


2 


(8) 

It is the Darcy’s law damping terms (as source term) that are 
added to the momentum equation due to phase change effect on 
convection. Where e =0.001 is a small computational constant 
used to avoid division by zero, the mushy zone constant C 


Pressure—Based Segregated Algorithm 


Update properties 


| 


Solve sequentially: 


U vel V vel W vel 


) 


Solve pressure—correction 


m~ 


(continuity) equation 


) 


Update mass flux, 


pressure, and velocity 


J 


Solve energy, species, 


turbulence, and other 


scalar equations 


describes the kinetic process in the mushy zone ordinary between 
10 and 10’. 


2.3. Fluent solver 
The Fluent software has two main solvers: the pressure-based 
solver and the density-based coupled solver (DBCS). Only the first 


method can be used to simulate the melting and solidification 


Pressure—Based Coupled Algorithm 


Update properties 


Solve simultaneously: 
equations and pressure-based 


continuity equation 


Update mass flux 


Solve energy, species, 


turbulence, and other 


scalar equations 


Fig. 2. Overview of the pressure-based solution methods [25]. 
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o Grid Independence 
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Fig. 3. Solution procedure overview for Fluent [30]. 
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problems. The pressure-based solver employs an algorithm that 
belongs to a general class of methods called the projection 
method [29] wherein the constraint of the mass conservation 
(continuity) of the velocity field is achieved by solving a pressure 
(or pressure correction) equation. The pressure equation is 
derived from the continuity and the momentum equations in 
such a way that the velocity field, corrected by the pressure, 
satisfies the continuity. The governing equations are nonlinear 
and coupled to one another; thus, the solution process involves 
iterations where the entire set of governing equations is solved 
repeatedly until the solution converges. Available in Fluent are 
two pressure-based solver algorithms, namely, a segregated 
algorithm and a coupled algorithm Fig. 2. 

Different discretization (interpolation methods) schemes 
are available for the convection terms in Fluent, including the 
first-order upwind scheme, power law scheme, second-order 
upwind scheme, central-differencing scheme, and the quadratic 
upwind interpolation scheme. The first-order upwind, power law, 
and second-order upwind schemes are mostly used with solidifi- 
cation and melting problems. The last two methods are more 
accurate than the first methods. In addition, the interpolation 
method of the face pressure for the momentum equations are the 
semi-implicit pressure-linked equation (SIMPLE), the SIMPLE- 
consistent (SIMPLEC), the pressure-implicit with splitting of opera- 
tors (PISO), and the fractional step method (FSM). Fig. 3 shows the 
solution procedure overview. More details about the solution, 
initialization, and discretization methods can be found in Ref. [25]. 

The physical properties of materials, such as density, 
thermal conductivity, heat capacity, and viscosity, may be tem- 
perature-dependent and/or composition-dependent. The temperature 
dependence is based on a polynomial, piecewise-linear, or 


a 


Melt fraction 


0O 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 
Time, min 


Melt fraction 


Time, min 


Fig. 4. Melt fraction vs. time for various values of the mushy zone constant with 


different models [31]. 


piecewise-polynomial function. Individual component properties 
are either defined by the user or computed via kinetic theory. 
Nevertheless, these physical properties can be defined as a 
constant value, a temperature-dependent function, or a user- 
defined function (UDF) that can be written in a specific program- 
ming language to define the temperature-dependence of the 
thermophysical properties. 

Some researchers deduced that the thermophysical properties 
of PCMs, such as density and viscosity, are dependent on tem- 
perature changes and are determined by specific correlations. 


p=p/(BT-T) +1) (9) 


p= exp(A+B/T) (10) 


Where p; is the density of PCM at the melting temperature T, 
f is the thermal expansion coefficient, and A, B are constant 
coefficients. The relationship between PCM and air is described by 
a specific model called the volume of fluid (VOF). This model 
defines the PCM-air system with a moving internal interface, but 
without inter-penetration of the two-phase fluid. Three condi- 
tions reflect the fluid’s volume fraction in the computational cells, 
denoted as &n: if %,=0, then the cell is empty of the nth fluid; if 
O&n=1, then the cell is full of the nth fluid; and if 0 <a, <1, then 
the cell contains the interface between the nth fluid and one or 
more other fluids. 

Shmueli et al. [31] numerically investigated the melting of 
PCM in a vertical cylinder tube with a diameter of 3 and 4cm, 
exposed to air above, and insulated at the bottom. The wall 
temperature was between 10 and 30°C above the melting 
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=+ Experimental 
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Fig. 5. Comparison between pressure discretization schemes and experimental 
results with different models [31]. 
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temperature of the PCM. They studied the effect of various 
numerical solution parameters such as the pressure velocity 
scheme (PISO and SIMPLE) and the pressure discretization regime 
(e.g., PRESTO and the weighted force method). They analyzed the 
influence of the mushy zone constant to the melting process, 
as shown in Fig. 4 and Fig. 5. They also reported no difference in 
the pressure velocity scheme (PISO and SIMPLE) as well as a 
considerable difference in the result for pressure discretization. 


3. Numerical simulation of PCM in thermal storage 
geometries 


3.1. Spherical capsule geometry 


LHTES in a spherical container represents an important case 
for thermal storage to be used in different engineering applica- 
tions such as in packed-bed storage. Moreover, most commercial 
companies produce the spherical capsule for this application. 

Assis et al. [32] presented a numerical and an experimental 
investigation of the melting process of RT27 with a PCM ina 
spherical geometry filled with 98.5% solid PCM as an initial 
condition of the simulation. One percent was left to account for 
the increase of the PCM volume during the phase change transi- 
tion. The variable density was defined during the liquid state with 
linear variation in the “mushy” state, and the VOF was used for 
the PCM-air system. The numerical model was based on the axial 
symmetry of the physical model using Fluent 6.0. Different design 
and operation parameters were studied, and the results indicated 
that the melting process of the PCM is affected by its thermal and 
geometrical parameters, including the Stefan number and the 
shell diameter. Assis et al. [33] numerically and experimentally 
studied the solidification process of RT27 as a PCM in a spherical 
shell filled with 98.5% liquid PCM. Different shell diameters were 
examined, and a transient numerical model was developed using 
Fluent 6.2; the numerical model was the same one used in the 
melting process of PCM [32]. 

Hosseinizadeh et al. [34] numerically studied the effect of 
various volume fractions of nanopractical copper on an uncon- 
strained melting rate in a spherical container. Different volumes 
of nanopractical copper (0, 0.02, 0.04) per volume were used, and 
85% of the PCM was filled. They used Fluent to develop an axi- 
symmetric numerical model and added the Darcy law to the 
momentum equation to account for the effect of phase change on 
convection. The VOF used for the PCM-air system and the power 
law scheme and the PISO method used for the pressure-velocity 
coupling were adopted to solve the momentum and energy 
equations, respectively. The PRESTO scheme was used for the 
pressure correction equation. Three different Stefan numbers 
were analyzed; the validation of the model was based on the 
experimental work by Assis et al. [32]. The results indicated that 
the nanopractical copper caused an increase in the thermal 
conductivity of the nano-enhanced PCMs compared with the 
conventional PCM. Unfortunately, the latent heat of fusions 
decreased. 

Tan et al. [35] reported an experimental and numerical work 
on the constrained melting of PCM inside a transparent spherical 
glass capsule. The Darcy law was added to the momentum 
equation to account for the effect of phase change on convection. 
They used Fluent 6.2.16 for the numerical simulation and the 
SIMPLE method for solving the governing equations. The power 
law differencing scheme was used to solve the momentum and 
energy equations, whereas the PRESTO scheme was adopted for 
the pressure correction equation. They reported that waviness 
and excessive melting of the bottom of the PCM was under- 
estimated by the experimental observation. This discrepancy 


was linked to the use of a support structure to hold the sphere 
that could have inhibited heat from reaching the bottom of the 
sphere. 

Xia et al. [36] developed an effective packed-bed thermal 
energy storage that contains a spherical PCM capsule. They 
studied the effect of the arrangement of the PCM spheres and 
the encapsulation of the PCM on the heat transfer performance of 
the LTES bed. They developed a 2D numerical model via Fluent 
6.2 using the standard k-€ model for turbulent flow; a SIMPLE 
algorithm was adopted for pressure and velocity coupling. They 
reported that random packing was more favorable than special 
packing for heat storage and retrieval; both the material and the 
thickness of the encapsulation had apparent effects on the heat 
transfer performance of the LTES bed. 


3.2. Square and rectangular geometry 


Arasu and Mujumdar [37] presented a numerical investigation 
on the melting of paraffin wax dispersed with aluminum (Al203) 
as a nano-PCM in a square enclosure heated from the bottom side 
and from a vertical side with different volume percentages 
(2% and 5%). They developed a numerical model using Fluent 
6.3.26. The computational domain was resolved with a fine mesh 
near the hot and cold wall to resolve the boundary layer; an 
increasingly coarser mesh was used in the rest of the domain to 
reduce computational time. They used a UDF to define the 
temperature-dependence of the thermophysical properties of 
PCM, wherein the first-order upwind difference scheme was used 
to solve the momentum and energy equations, and the PRESTO 
scheme was adopted for the pressure correction equation. They 
reported that the effective thermal conductivity of a paraffin wax 
latent heat storage medium could be increased significantly 
by smaller volumetric concentrations of alumina particles in the 
paraffin wax. 

Pinelli and Piva [38] developed a 2D numerical model using 
Fluent to study the effects of natural convection on the PCM 
process in terms of temperature distributions, interface displace- 
ment, and energy storage, and to enhance the agreement between 
the experimental work with the numerical one when the con- 
duction dominated is considered only for the heat transfer of the 
PCM. The mushy zone constant for the model was approximately 
1.6*10°, and the cylinder cavity was filled with n-octadecane as 
PCM in the solid state and heated from above with an electrical 
heater. The overall heat transfer coefficient U considered conduc- 
tion through the polystyrene layer as well as the convection and 
radiation heat exchanged with the environment. The model was 
validated with an analytical solution wherein the conduction- 
dominated freezing and numerical methods were considered in 
the convection-conduction. They reported that the agreement 
between numerical and experimental results significantly improved 
when the presence of circulation in the liquid phase, instead of in 
the conduction-dominated process only, was considered. 

Researchers reported that the heat transfer in PCM could 
be enhanced significantly by adding nanopractical copper to the 
PCM (nanofluid). Such addition improves the PCM’s thermal 
conductivity. Wu et al. [39] numerically investigated the melting 
processes of Cu/paraffin nanofluid PCMs. They developed a 2D 
numerical model using Fluent 6.2. The enclosed cavity was filled 
with 2.5-cm high PCM and 1 cm air; the top of the cavity was 
enclosed and taken by air. The VOF was adopted to describe the 
relationship between the PCM and air with a moving interface, 
but without inter-penetration of the two media. The QUICK 
differencing scheme was used to solve the momentum and energy 
equations, whereas the PRESTO scheme was adopted for the 
pressure correction equation. They reported a 13.1% decrease in 
the paraffin’s melting time. 
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Li et al. [40] numerically studied the melting process of a PCM 
in concentric horizontal annuli of a square external shell. A 2D 
numerical model was adopted using Fluent 6.0. They reported 
that the top part of the PCM melts faster than the bottom part as a 
result of the convection-dominated mode, which accelerates the 
melting of the top part of the annulus. 

Khodadadi and Hosseinizadeh [41] reported that the disper- 
sion of nanoparticles in PCMs enhanced the material’s thermal 
conductivity in contrast to the base material. They used a 2D 
numerical model within the commercial code of Fluent 6.2.1. The 
SIMPLE algorithm and the QUICK differencing scheme were used 
to solve the momentum and energy equations, whereas the 
PRESTO scheme was adopted for the pressure correction equation. 
The Darcy Law damping terms were added to the momentum 
equation. The results indicated an increase in the thermal con- 
ductivity of the nano-enhanced PCMs in contrast to the conven- 
tional PCM. Unfortunately, the latent heat of fusion decreased. 


3.3. Cylindrical capsule geometry 


Cylindrical geometries are considered most promising for 
devices for commercial heat exchangers, such as the double pipe 
heat exchanger and shell and the tube heat exchanger, because of 
their high-efficiency in a minimum volume. Most researchers use 
this geometry for LHTES, filling the PCM in the tube side or in the 
annulus (shell). Commercial products of PCMs are delivered in 
these geometries. 

Darzi et al. [42] used Fluent 6.2 to simulate the melting 
process of the n-octadecane as a PCM in the concentric and 
eccentric double pipe heat exchanger. They studied the effect of 
the internal tube position on the melting rate. They developed a 
2D numerical model (R, 0) and adapted low power to solve the 
momentum and energy equations. The SIMPLE method was used 
for pressure-velocity coupling, and the PRESTO scheme was used 
for the pressure correction equation. The results indicated that 
the melting rate was the same in the early stage of melting, and 
then decreased in the concentric model. 

Gou and Zhang [43] reported that the solidification process 
of KNO3-NaNO3 as PCM in solar power generation using the 
direct steam technology is controlled by pure conduction. They 
neglected the influences of natural convection when comparing it 
with heat conduction. They developed a 2D numerical model via 
Fluent 6.2 and studied the storage with and without foil based on 
the same tube diameters and PCM. The results indicated how the 
heat transfer and discharge time were affected by the changes in 
the geometry of the aluminum foil and in the tube radius. They 
reported that the discharge process significantly accelerated with 
the addition of the aluminum foil. 

Colella et al. [2] designed a medium-scale LHTES unit for 
district heating systems. The shell-and-tube LHTES that included 
RT100 as a PCM in the shell side was used to transfer heat from 
the district to the heating networks of the building; the heat 
transfer fluid was water. The thermal storage was charged during 
the night using CHP plants and then discharged during the day at 
the peak of the thermal request. A 2D numerical model (R, X) was 
developed using the Fluent simulation software to study the 
thermal behavior of the PCM and the HTF. Grid meshing of the 
computational model was refined systematically. The SIMPLE 
algorithm was adopted to solve the pressure-velocity coupling. 
The second-order upwind scheme was used for the convective 
fluxes. A naturally expanded graphite matrix with 15% volume 
was added to the paraffin to enhance the thermal conductivity of 
the composite (paraffin-graphite) to achieve the heat fluxes 
required for the medium scale application. Different operation 
strategies and parameters were studied, including the time-wise 
variation of the liquid fraction, PCM temperature, and the HTF 


outlet temperature. Based on their analysis, the best options for 
energy storage density were characterized by latent storage units 
installed on the heating network of the building rather than on 
the primary district heating network. 

Buyruk et al. [44] studied the solidification process of water 
around different staggered and inline cylinders placed in a 
rectangular ice storage tank. They developed a 2D numerical 
model using Fluent to depict the temperature distributions and 
ice formation in the tank. They adopted Darcy’s law and the 
Kozeny-—Carman equation to model the flow and permeability of 
the mushy zone, respectively. Based on a fixed ice storage tank 
volume, the solidified area of four inline cylinders was larger than 
that for six inline cylinders. 


4. Numerical simulation of the PCM applications 
4.1. Numerical simulation of PCM for electronic devices applications 


Shatikan et al. [45] numerically presented the melting process 
of a PCM with internal fins exposed to air from the top and heated 
from the horizontal base. They investigated different fin dimen- 
sions with a constant ratio between the fin and the PCM thick- 
ness. They used Fluent 6.0 to develop 2D and 3D numerical 
models and adopted the VOF model to describe the relationship 
between the air and the PCM; the SIMPLE algorithm was used for 
the pressure-velocity coupling. The 2D model was considered 
because the 3D model was time consuming. They reported that 
the results for the 2D and 3D models were equal in some 
dimensions, and the transient phase change process of the PCM 
depended on the operation and design parameters of the system. 
Such parameters were related to the temperature difference 
between the base and the mean melting temperature, as well as 
to the thickness and height of the fins. Shatikian et al. [10] 
analyzed the same model with constant heat flux applied from 
the horizontal base. 

Wang and Yang [11] numerically studied the cooling technol- 
ogy of portable handheld electronic devices using PCM. They 
developed half computational grids with the 3D numerical model 
using Fluent 12 to investigate the effect of a different amount of 
fins and various power’s heating level. They adopted the VOF 
model to describe the relationship between the air and the PCM. 
In the VOF model, the Geo-Reconstruct, which is the volume 
fraction discretization scheme, was used to calculate the transient 
moving boundary. If the transient temperature of the calculation 
cell in the PCM domain was equal to or higher than the melting 
temperature, the continuity and momentum equations were 
calculated in the system; the SIMPLE algorithm was used for the 
pressure-velocity coupling. They studied the orientation of PCM 
with the heat sink; the model was validated using the method by 
Fok et al. [46]. The numerical results indicated that the orienta- 
tion had insufficient influence on the transient thermal perfor- 
mance of the hybrid cooling system, and that the transient surface 
temperature caused a discrepancy of around 8.3%-10.7%. They 
reported that PCM in aluminum heat sink with six fins could 
provide stable temperature control and good cooling ability. 

Yang and Wang [12] numerically investigated the cooling 
technique using a hybrid PCM-based heat sink within a closed 
system to absorb constant and uniform volumetric heat generated 
from portable handheld electronic devices. They tested different 
computation grids and time steps to validate the model with 
the experimental work; various power levels (2-4 W), different 
orientations (vertical, horizontal, slanted), as well as charging 
and discharging times were conducted in the simulation model. 
The melting rate increased when a higher power level was 
transferred to the system. They reported that the maximum 
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temperature of the hybrid system using PCM could be well 
controlled under 320 K if the melting ratio was under 0.9 or if 
the melting time was shorter than 6801.2 s. 

Ye et al. [47] figured out the correlation between influential 
variables, such as the liquid fraction and thermal storage time, the 
transient heat flux and melting period, and the solid fraction and 
the solidification time by numerically studying the heat transfer 
and fluid flow in a plate fin filled with PCM for rapid heat 
absorbed/released techniques. They developed a 2D numerical 
model through Fluent. The cavity was filled with 85% PCM and 
15% air to consider thermal expansion during the solid-liquid 
phase change; finely structured meshing grids were installed near 
the heating and cooling plate wall to resolve the boundary layers. 
The VOF model was adopted to describe the relationship between 
the air and the PCM. The UDF was written in C language to 
account for the temperature dependence of the thermophysical 
properties. The SIMPLE algorithm was considered for the pressure- 
velocity coupling, and different time steps were used for the 
energy storage and energy release. Exploring via a 3D numerical 
model, they reported no difference in the results between the 2D 
and 3D models. The correlations they obtained will be useful for 
future component design and system optimization. 

Ye et al. [48] performed further studies on the effect of PCM 
cavity volume fraction on the heat transfer behavior and fluid 
flow in LHTES. Different parameters of LHTES performance, such 
as the volume expansion ratio, the time of complete thermal 
storage, heat flux, liquid fraction, and velocity and temperature 
fields, were investigated. They used the same model in Ref. [47], 
except that the time step was chosen in relation to the grid size 
via the Courant-Friedriches-Lewy number. They reported that 
the volume expansion ratio decreased as cavity volume fraction 
x increased. By contrast, the time for complete energy storage 
increased. 

Kandasamy et al. [49] numerically and experimentally inves- 
tigated the use of a PCM-based heat sink for quad flat package 
(QAP) electronic devices. They developed one quarter of the 
physical model as well as a 3D numerical model using the 
commercial Fluent 6.2. The VOF was used to describe the relation- 
ship between the air and the PCM. The PISO scheme was used for 
the pressure-velocity coupling. The cooling performance of the 
PCM-based heat sinks improved, compared with the case without 
PCM. The full melting of the PCM and the melting rate also 
increased as the input power increased. 

Sabbah et al. [50] presented the thermal and hydraulic 
behavior of a micro-channel heat exchanger with MCPCM water 
slurry for heat dissipation of a high-power electronic device. They 
developed a 3D numerical model using Fluent 6.2, and the SIMPLE 
algorithm was used for velocity-pressure coupling. The heat 
transfer in both the cooling fluid and the metal heat sink on the 
overall performance of the system was considered. The results 
showed a significant increase in the heat transfer coefficient that 
depended on the channel inlet and outlet and on the selected 
MCPCM melting temperatures. Lower and more uniform tem- 
peratures throughout the electronic device could be achieved 
using less pumping power compared with the case using only 
water as the cooling fluid. 

Hosseinizadeh et al. [51] numerically and experimentally 
presented the thermal management of electronic products using 
PCM. They compared the effect of the PCM cooling technology in 
contrast to the normal ones. Different design parameters were 
studied to determine their effects on the temperature of the heat 
sink and on the melting of PCM inside the heat sink. 85% of the 
heat sink height was filled with PCM, and 15% of the height was 
modeled in the air region for the solid-to-liquid conversion in the 
PCM expansion. A 2D numerical model was adopted via Fluent. 
Density was considered constant in the solid state, whereas liquid 


density was a function of the temperature and the thermal 
expansion coefficient. The VOF was adopted to describe the 
PCM-air system. The PISO algorithm was used for the pressure- 
velocity coupling. They set the time step at 0.0001 s at the early 
stage of simulation and changed it to 0.001s after 3 min of 
simulation. The increased fin thickness showed only a slight 
improvement in thermal performance, whereas an increase in 
the number of fins and fin height resulted in an appreciable 
increase in the overall thermal performance. 


4.2. Numerical simulation of PCM for HVAC equipment 


Tay et al. [52] modeled different tube configurations in a shell 
and tube heat exchanger with PCM in the shell side. They 
developed a 3D CFD model using ANSYS code to analyze the 
transient heat transfer during the phase change process. Three 
grid resolutions were evaluated for the tube, HTF, and PCM. 
ANSYS CFX 12.1 was used to complicate the physical properties 
of the mixture. Three domains were created using CFX.PRE 
version 12.1, namely, liquid for HTF, liquid for PCM, and solid 
for the tube. An experimental apparatus was developed to 
compare the predicted simulation results of both freezing and 
melting processes. They found that the results of the CFD melting 
model were generally longer than the experimental results. The 
thermal behavior of the freezing and melting processes of PCM 
was different from those in the experimental results when the 
natural convection was ignored. 

Antony and Velraj [53] studied the heat transfer and fluid flow 
through a tube in a regenerative PCM heat exchanger (shell and 
tube) module incorporated into a ventilation system. The PCM 
was encapsulated in the shell side of the module for a free cooling 
system. Fluent software was used to analyze the heat transfer and 
flow through the model, which consisted of an air flow passenger 
and two air spacers on top and at the bottom. PISO algorithm was 
adopted for the pressure-velocity coupling. The second-order 
upwind scheme was used for the approximation of the convective 
terms. The k-€ model was used for the turbulence modeling of 
the heat transfer fluid. This module was validated experimentally. 
The results from the experiment and the calculated results were 
compared. 

Xia and Zhang [54] prepared an acetamide/expanded graphite 
composite (AC/EG) with 10% (mass fraction) as PCM to be used in 
active solar applications, such as solar-driven solid/liquid desic- 
cant dehumidification systems and solar-driven adsorption/ 
absorption refrigeration systems that require a higher heat source 
temperature of over 60 °C. They measured the thermal properties, 
including melting/freezing points, thermal conductivities, and the 
heat of fusion, using the transient hot-wire method and a 
differential scanning calorimeter. Two-thirds of the test model 
was filled with solid PCM for the thermal expansion during the 
solid-liquid phase change. They developed a 2D numerical model 
using Fluent 6.2 to study the release and absorption of heat 
in LHTES. The numerical model was validated experimentally. 
The thermal conductivity of the composite AC/EG (2.61) was six 
times that of pure metal; 0.043 for solid at 60 °C and 0.2 for liquid 
at 90°C. The phase change temperature shifted and the latent 
heat decreased compared with that of pure AC. The heat storage 
and retrieval durations of the LTES unit with the AC/EG composite 
showed 45% and 78% reductions, respectively, compared with 
those for pure AC. 

New internal fins inside spherical and cylindrical PCM encap- 
sulation in the HVAC system for buildings were numerically 
investigated by Siva et al. [55]. They developed a numerical 
model using Fluent 6.3 to study the charging and discharging 
process of PCMs for the two configurations. Incorporating fins for 
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different configurations reduced the total solidification time by 
65%-72%. 


4.3. Numerical simulation of PCM for building applications 


LHTES is used extensively in building structures to remove 
thermal load and to improve the thermal comfort of occupants. 
PCM can be found in different items of the building structure, 
such as the walls, roof, floor, ceiling, window shutters, plaster- 
board, and tiles. Most researchers mentioned in the present study 
used CFD to simulate the process, reporting a good agreement 
between the simulation and the experimental results. 

Silva et al. [56] incorporated a PCM into a typical clay brick 
masonry enclosure wall in Portugal as a passive technology to 
enhance the building’s energy efficiency. They experimentally and 
numerically evaluated how the PCM reduced internal tempera- 
ture fluctuations and increased the time delay between internal 
and external conditions. They generated the mesh grid in CFD and 
developed a 2D numerical model in ANSYS Fluent 12. The SIMPLE 
algorithm was used to solve the governing equations. The second- 
order upwind scheme was used to solve the pressure, momen- 
tum, and energy equations. They reported that incorporating the 
PCM into the wall contributed to the attenuation of the indoor 
space temperature swing, which reduced from 10 to 5 °C, and to 
the thermal amplitude. The time delay also increased by approxi- 
mately 3 h. 

Joulin et al. [6] analyzed the energy conservation of a building 
with PCM conditioned in a parallel epipedic polyfin envelope to 
be used in passive solar walls. They developed a 1D Fortran 
numerical model and a 2D model using commercial Fluent to 
study the thermal behavior of the PCM to validate the experi- 
mental work. The self-developed numerical model was based on 
the enthalpy formulation and on the fully implicit finite differ- 
ence solution method used for discretization. The equations were 
solved using TDMA. The liquid fraction was updated in each 
iteration until convergence was achieved. For the commercial 
Fluent, the SIMPLE method was used to solve the governing 
equations. The power law differencing scheme was used to solve 
the momentum and energy equations, whereas the PRESTO 
scheme was adopted for the pressure correction equation. They 
reported that the 1D and 2D numerical simulations did not satisfy 
because the results were very far from the experimental ones: the 
codes did not represent the supercooling phenomenon. Numerical 
modeling showed that the supercooling phenomenon must be 
taken into account correctly to predict the PCM thermal behavior. 

Susman et al. [57] constructed four 150 mm? prototype PCM 
sails from paraffin/low-density-polyethylene composites installed 
below the ceiling of an occupied office space in London in the 
summer. During daytime, the modules were hung on an internal 
rig 2.5m above the ground; after 8 p.m., the modules were 
transferred to the outside rig to discharge the energy, and were 
moved to the inside rig at 7 a.m. They developed a semi-empirical 
model of the modules in Fluent using an enthalpy-porosity 
formulation to model phase change. The modules were validated 
well using the temperature measurements, including a notable 
divergence when the maximum liquid fraction was reached. 


4.4. Numerical simulation for other applications 


Bo et al. [58] introduced a metal PCM energy storage in a 
high-temperature heat pipe steam boiler used for direct steam 
generation in a solar thermal power generation system. The Al-Si 
eutectic alloy contained 12.07% Si used as a PCM at a melting 
temperature of 577 °C The temperature distributions of the heat 
reservoir tank and the heat storage medium under high solar heat 
flux were analyzed using a 2D numerical model via Fluent. The 


newly developed heat storage boiler could withstand a focused 
solar energy flux of 400 kW/m? and was compatible with the heat 
storage medium. 

Tan et al. [16] used water as a PCM to recover and store cold 
energy from a cryogenic gas in a cryogenic cold energy system. 
They developed a 2D numerical model using Fluent to analyze the 
effect of the thermal boundary and thermal resistance on the 
freezing process in an LHTES system. Their numerical model was 
based on half the domain because the geometrical model was 
considered symmetric. The standard k-€ model was adopted to 
calculate the internal forced convection of the gaseous coolant. 
The CFD model was validated experimentally and showed good 
agreement with the experimental results. They reported that the 
dimensionless numbers, such as the Biot number and the Stefan 
number of the PCM, and the Stanton number of the coolant 
flowing in the tube had a noticeable influence on the PCM 
freezing characteristics. A larger Biot number was beneficial to 
the heat transfer between the PCM and the coolant, promoting a 
higher freezing rate. A higher Stefan number appeared to result in 
larger freezing rates in a fixed axial position. Moreover, the frozen 
layer slope along the tube was steeper at larger Stanton numbers. 

Xiaohong et al. [59] numerically studied a PCM at high melting 
temperatures (1120 K) in a heat pipe receiver in an advanced 
solar dynamic power generation system. A 2D numerical model 
was adopted via Fluent 6.2 to simulate the heat transfer of the 
PCM canister of the heat pipe receiver. UDF was used in the 
simulation. The numerical results were compared with NASA’s 
numerical results. The temperature of the heat pipe wall and 
the outer wall of the PCM canister in the numerical model and in 
the NASA simulation differed only slightly. These results can be 
used to guide the design and optimization of PCM canisters for 
heat pipe receivers. 

Koizumi and Jin [60] proposed a new compact slab container with 
an arc outer configuration to promote direct contact melting. They 
studied two configurations of the close contact melting process, 
namely, flat slab container and curved slab container. The 2D 
numerical model was performed via Fluent 6.3. The VOF model was 
adopted to describe the PCM-air system. A first-order accurate 
upwind difference scheme was adopted for the convection terms 
because the convection velocity in the liquid PCM was extremely 
small (within several mm/s). The simulated results quantitatively 
elucidated the experimental melting process from beginning to end. 

MacPhee and Dincer [61] numerically investigated and ther- 
modynamically analyzed the melting and solidification process of 
a de-ionized water as PCM in a spherical capsule geometry. They 
developed a 3D model using Fluent 6.0 and studied the effect of 
the HTF inlet temperature and the flow rate. The numerical model 
was validated and was found to be in good agreement with the 
experimental data. The variations of the incoming HTF tempera- 
ture had a much larger effect on the charging times and on the 
energy and exergy efficiencies compared with the changing of the 
HTF flow rate. 

Reddy [62] presented a solar-integrated collector storage 
water heating system with PCM to optimize the solar gain and 
heat loss characteristics through a 2D numerical model using 
commercial Fluent and UDF to express the variable heat flux 
condition on the absorber plate. Different numbers of fins were 
attached to the PCM water solar system to enhance the heat 
transfer. The latent heat storage with nine fins was optimal at 
maximum water temperature and had minimum heat loss to the 
surroundings. 

Rao et al. [63] discussed the use of PCM with an aging 
commercial LiFePO, power battery to study the thermal behavior 
of the PCM and the cell for thermal energy management. 
They developed a 3D computation model for a single cell and a 
battery package using the commercial computational Fluent. They 
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reported that the design for battery thermal management should 
include a proper selection of the thermal conductivities for the 
PCM and the cell because the thermal conductivity increased, 
whereas the latent heat of fusion decreased. 

Darkwa and Su [64] studied the effect of particle distribution on 
the thermal performance of micro-encapsulated PCMs in different 
configurations of a composite high conductivity laminated MCPCM 
board. They developed a 3D numerical model using Fluent 6.3. The 
SIMPLE algorithm was used for the pressure-velocity coupling, and 
the second-order upwind scheme was adopted for both momentum 
and energy equations. The PRESTO pressure discretization scheme 
was used for the pressure correction equation. The thermal response 
times for the rectangular and triangular geometries were approxi- 
mately half of that for the pyramidal geometry during the cooling 
and heating processes of the board. 

MacPhee et al. [65] numerically studied the effect of different 
capsule geometries, HTF inlet temperatures, and HTF flow rates 
on energy and exergy efficiencies for the solidification process in 
encapsulated ice thermal energy storage (EITES). They used Fluent 
6.3 to simulate different geometries (slab, cylindrical, spherical), 
three HFT flow rates, and three HTF inlet temperatures to achieve 
the highest efficient charge method for the thermal storage. They 
recommended the exergy efficiency methods to study system 
performance, stating that these methods provide better insight 
into system losses. They also advised EITES designers to increase 
both the flow rate and inlet HTF temperature to achieve full 
system charging in an acceptable amount of time. 


5. Conclusion 


The sustainable thermal energy storages for different engineering 
applications, such as building structures and HVAC equipment, are 
indispensable to reducing greenhouse-gas emission. LHTES is an 
important component of thermal solar engineering systems, playing 
a key role in improving the energy efficiency of these systems. 
Thermal energy storage systems, especially LHTES, have gained 
widespread attention in relation to global environmental problems 
and energy-efficiency improvement. The following conclusions can be 
summarized for CFD applications in latent heat thermal storage. 


1. The numerical solution for the PCM phenomena is more 
accurate than the analytical solution and can be used in 
different engineering conditions. 

2. The CFD software is used successively to simulate the 
application of PCM in different engineering applications, 
such as electronic cooling technology, building thermal 
storages, and HVAC. 

3. Different CFD software are used in PCM engineering, including 
ANSYS Fluent, Comsol Multiphysics, and Star-CCM+. The 
most commonly used software is ANSYS Fluent. 

4. The numerical results of the 2D numerical model are generally 
the same as those of the 3D numerical model. Researchers 
reported that the time required for simulation was reduced. 

5. The application of CFD in designing PCM thermal storages 
is a feasible method because of the highly accurate results 
of the experimental work. CFD also delivers optimization 
tools to help users achieve maximum efficiency while 
saving time and money. 
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